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Abstract 

Suppose a coin with unknown probability p of heads can be flipped 
as often as desired. A Bernoulli factory for a function / is an algo¬ 
rithm that uses flips of the coin together with auxiliary randomness to 
flip a single coin with probability f{p) of heads. Applications include 
perfect sampling from the stationary distribution of certain regener¬ 
ative processes. When / is analytic, the problem can be reduced to 
a Bernoulli factory of the form f{p) = Cp for constant C. Presented 
here is a new algorithm that for small values of Cp, requires roughly 
only C coin flips. From information theoretic considerations, this is 
also conjectured to be (to first order) the minimum number of flips 
needed by any such algorithm. 

For large values of Cp, the new algorithm can also be used to build 
a new Bernoulli factory that uses only 80% of the expected coin flips 
of the older method. In addition, the new method also applies to the 
more general problem of a linear multivariate Bernoulli factory, where 
there are k coins, the fcth coin has unknown probability pk of heads, 
and the goal is to simulate a coin flip with probability Cipi + - ■ ■+CkPk 
of heads. 
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1 Introduction 


The n otion of a Bernoulli factory was introduced in lAsmussen. Glvnn. a.nd Thorisson 
(119921) in the context of generating samples exactly from the stationary dis¬ 
tribution of a regenerative Markov process. A Bernoulli factory works as 
follows. Suppose we have the ability to draw independent identically dis¬ 
tributed (iid) Bernoulli random variables, each of which is 1 with probability 
p and 0 with probability 1 —p (write X ~ Bern(p).) Then given a function 
/, the goal is to use a random number of draws from X to build a new ran¬ 
dom variable which is also Bernoulli, but with chance fjp) of b eing 1 for a 
specihed function /. In Asmussen. Glvnn. a.nd ThorissonI f 1992[l . the needed 
function was a linear function, namel y a constant times p. This simple case 
generalizes: in iNacu and Peres! (1200511 it was shown that the ability to draw 
from /(p) = 2p could be used to build a Bernoulli factory for any analytic / 
that was bounded away from 1. 

The focus here is on building a nearly optimal linear Bernoulli factory 
where Cp is known to be small. This is nearly optimal in the sense that 
it uses (to hrst order) only C flips of the coin, and th ere is s t rong evide nce 
to indicate that at least C flips are necessary. As in iHuberl (Ito appearl) . a 
Bernoulli factory can be dehned as follows. 


Definition 1. Given p* e (0,1] and a function / : [0,p*] —)■ [0,1], let A 
be a computable function that takes as input a number u G [0,1] together 
with a sequence of values in {0,1}, and returns an output in {0,1}. For any 
p G [0,p*], Xi,X 2 ,... iid Bern(p), and U ~ Unif([0, 1]), let T be the infimum 
of times t such that the value of A{U, Xi, X 2 ,...) only depends on the values 
of Xi,..., Xt- If the following holds, then call A a Bernoulli factory. 


1. T is a stopping time with respect to the natural filtration that is finite 
with probability 1. 


2. A([/,Xi,X 2 ,...)~Bern(/(p)). 

Call T the running time of the Bernoulli factory. 

Colloquially, a draw X ~ Bern(p) will be refereed to as a coin flip, or 
more specihcally, a p-coin flip. The result X = 1 corresponds to heads on 
the coin, while X = 0 indicates tails. So a Bernoulli factory attempts to flip 
a coin with /(p) chance of heads, by using a random number of coin flips 
from the original coin together with some auxiliary randomness. 
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Asmussen, Glynn, and Thorissoru fjl992f l introduced Bernoulli factories 
for an application in perf ect simulation, but did not show that they exist. 
Keane and O’Brien f 1994|) constructed the hrst general Bernoulli factories, 


showing that such a factory with hnite running time existed if and only if 
f{p) was continuous over [0,p*] for some p* G (0,1], and either it holds that 
f{p) is identically 0 or 1, or that both /(p) and 1 — /(p) are polynomially 
bounded away from 0 and 1 over the allowable range of p. 

Their strategy for building a Bernoulli factory was to construct Bernstein 
polynomials that approximated the function /(p) as closely as possible. Bern¬ 
stein polynomials are linear combinations of functions of the form p^(l— 
where both n and k < n are nonnegative integers. Such polynomials can be 
created from the coin by flipping it n times and seeing if exactly k heads 
and n — k tails result. Keane and O’Brien could show that the running time 
T was hnite with probability 1 for their algorithm, but not much more. In 
particular, they could not show any bounds on the average running time, or 
eve n that it was hii i te. 

Nacu and PeresI ( 2005 1 developed this approach further, and showed that 


Bernstein polynomials could be constructed tightly enough that the running 
time would have a hnite expectation. In addition, they showed that the tail 
of the distribution of their running time declined exponentially. Moreover, 
their work contained a proof that /(p) = 2p is in a sense the most important 
function, since it can be used to construct a Bernoulli factory for any function 
that is both real analytic over [0,1] and bounded away from 1. 

However, their approach was not a practical algorithm. While the num¬ 
ber of coin hips had hnite expectation, the amount of memory and time 
needed to com pute the function A grew exponentially with the number o f 
hips. Work of iLatuszvhski. Kosmidis, Papspiliopoulos. and Roberts! (1201111 
solved this issue, and gave the hrst practical implementation of the Nacu 
and Peres approach. Their approach created a pair of reverse time processes, 
one a supermartingale, the other a submartingale, that converged on the tar¬ 
get /(p). The values could be computed without the exponential overhead 
associated with the Nacu-Peres algorithm. 

To bound /(p) = Cp away from 1, they considered the function /(p) = 
min{Cp, 1 — e} so that the function was dehned over the entirety of [0,1]. 
However, this was not strictly necessa ry, as in the original application of 


Asmussen. Glvnn. and Thorissoru (1199211 , it was possible to easily insure that 


f{p) < 1—e. By not trying to sample from the function /(p) = min{Cp, 1—e} 
for all values of p, but only for those with Cp < 1 — e, a new approach became 
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possible. 

This new approach was developed in lHuberl fjto appeaii ). and was the first 
that did not begin with Bernstein polynomial approximations. Instead, this 
approach used flips of the coin to alter the problem in ways that insured that 
the hnal output had the correct distribution. For instance, suppose the goal 
was to generate a coin with probability of heads 2p. Then flip the original 
coin once. If the coin is heads, the the output is heads. Otherwise, it is 
necessary to flip a p/(l — p)-coin. 

That way, the chance the final output is heads is p(l) + (1 —p)-pj (1 —p) = 
2p. By advancing carefully in this manner, it was shown how to build a 
Bernoulli factory such that for Cp < 1 — e where e is a known constant. 


E[T] < g.SOe 


-1 


Moreover, the same work showe d that this running time is the best possi¬ 
ble up to a constant. Specihcally, in lHuberl fjto appearl ) it was shown that that 
any Bernoulli factory (to hrst order) must use on average at least 0.0406“^ 
coin flips. It remains an open question what the best constants for the 
lower and upper bounds are, although there is strong reason to believe (see 
Section 0]) that on average at least C flips (to first order) are necessary to 
generate a Cp coin. 

The primary application of the Bernoulli factory, starting with Asmussen, 
Glynn, and Thorisson (1992), is to generate perfe ct samples from the station¬ 
ary di stribution of regenerative Markov chains. In Lee. Dou c et and Latusvznskil 
( 2014 1. it was shown how to use the Bernoulli factory in Huberl fjto appear l 
to generate coins where Cp <2/3. Under this condition, the older algorithm 
gave a bound of 38 coin flips on the average number needed. 

Without going into the details of their algorithm, by tripling the expected 
running time of there algorithm, it is possible to ensure that Cp < 2/9. 
Under these conditions, the expected number of flips for the new algorithm 
is bounded above by 7.8 (see Theorem [1]) giving an algorithm that only 
requires 62% as many work on average as the old one. 

This work gives the following results. 


1. For Cp small, an algorithm will be given that uses only C coin flips on 
average. 

2. For Cp at most 1 — e for known e, an algorithm will be given that uses 
only 7.57C'e“^ coin flips on average. 
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3. The new algorithm can be extended from the function Cp for single 
variate coins to the multivariate coin problem where there are k coins 
with unknown means pi,..., pk- Suppose the goal is to generate Bern(r) 
where 

r{pi,...,pk)=Cipi^ - hCkPk- (1) 

Then setting C = Ci + ■ ■ ■ + Ck, the algorithm in the multivariate case 
has running time equal to the single coin case. 

More precisely, the running time of the new algorithm is given as follows. 

Theorem 1. Suppose it is known that Cp < M for a constant M < 1/2. 
Then there exists an algorithm for producing a Cp-coin that uses on average 
at most 

C , 

(1-2M)(1 + C'p) ^ ^ 

coin flips. 

This theorem is shown in Section 13.51 The multivariate version is similar, 
and is shown in Section [5l 

Theorem 2. Let C = Ci + ■ ■ ■ + Ck, and r = Cipi + ■ • • + CkPk, Suppose it is 
known that r < M for a constant M < 1/2. Then there exists an algorithm 
for producing an r-coin that uses on average at most 

C , 

(l-2M)(l + r) ^ 

flips from among the k coins. 

The remainder of this paper is organized as follows. Section [2]presents the 
algorithm for small r, and shows correctness and the bound on the running 
time. Section inigives the extension to the multivariate problem and for larger 
values of r, and also includes the proofs of correctness and the bound on the 
running time. Finally, Section 0] considers why C flips is likely the best 
possible. 


C- 


15.2 


1 — 2M + r 


C 


15.2 


1-2M + Cp 
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2 The algorithm for small Cp 


Let r = Cp. The first piece of the algorithm is a method for drawing from 
the logistic Bernoulli factory 


that uses T coins, where E[T] = C/{I + r). 

As usual, say that X is exponential with rate A (write X ~ Exp(A)) if X 
has density fx{,s) = exp(—As)l(s > 0). Here !(■) is the indicator function 
that evaluates to 1 when the argument is true and 0 when the argument is 
false. The following basic facts about exponentials will prove useful. 

Fact 1. Let X ~ Exp(Ai) and Y ~ Exp(A 2 ) he independent. Then P(X < 
y) = Ai/(Ai + A 2 ). 

Fact 2 (Memoryless). If X Exp{X), then for s > 0, the conditional dis¬ 
tribution of X — s given X > s is exponential with rate A as well. That is, 
[X -s|X > s] ~ Exp(A). 

Exponentials can be employed to dehne a one dimensional Poisson point 
process. 

Definition 2. Let Ai, A 2 ,.. .he independent and identically distributed (iid) 
exponential random variables with rate A. Then 

P = {Ai, Ai + A 2 , Ai + A 2 + A 3 ,...} 

forms a Poisson point process on [0, cxd) of rate A. For [a, b] C [0, cx)), Pn[a, b] 
is a Poisson point process on [a, b] of rate A. 

Several well known facts about Poisson point processes are useful. 

Fact 3. The converse of the definition holds: any Poisson point process 
P C [0, cxd) of rate A with points 0 < Pi < P 2 <■■ ■ has Pi ~ Exp(A) and 
Pi — Pi-i ~ Exp{\), and all these exponentials are independent. 

Fact 4. Let P = {Pi,P 2 ,---} be a Poisson point process. Let Bi,B 2 ,... 
he a sequence of iid Bern{p) random variables. Then P' = {P* : P* = 1} 
is a Poisson point process of rate \p. [The process P' is called the thinned 
process.] 
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Fact 5. Let Pi and P 2 be independent Poisson point processes of rate Ai and 
X 2 over [0, 00). Then Pi U P 2 is a Poisson point process of rate Ai + A 2 over 
[0, cx)). 

Fact 6. The expected number of points in a Poisson point process of rate A 
over [a, b] is Poisson distributed with mean \{b — a). 

These ideas can be used to build the logistic Bernoulli factory for r/(l+r). 


Logistic_Bernoulli_Factory Input: C 

1) X ^ 0, draw A <- Exp(l) 

2) Draw T Exp(C') 

3) While X = 0 and T < AI 

4) Draw B Bern(p) 

5) If i? = 1 then X = 1, else T <— T -I Exp(C') 

6) Return X 


Note that line 4 can be accomplished i n constant tim e (with 0(k) pre¬ 
processing time) using the Alias method of IWalkerl fjl974|) . 


Lemma 1. The output o/Logistic_Bernoulli_Factory is a Bernoulli with 
mean r/(l -F r). 


Proof. Let Ti, T 2 ,... be the successive values of T taken on in the algorithm, 
and Bi, B 2 ,... the successive values of B. Since Tj+i — Tj is Exp(C), the {Tj} 
form a Poisson point process P of rate C. Let P' be the points Ti & P with 
Bi = 1. Then P' is a point process with rate Cp = r. 

Let T[ = min{P'}. The while loop examines the P' process, and returns 
1 if T[ < A, and 0 otherwise. By Fact [3], T[ ~ Exp(r), and A ~ Exp(l), so 
P(T( < A) = r/(l -F r) by Fact [H □ 


Lemma 2. In one call to Logistic_Bernoulli_Factory, the expected num¬ 
ber of coin flips needed is C/{1 -F r). 

Proof. The Poisson process of rate C combined with the process of rate 1 
forms a Poisson point process of rate C + 1. The chance that any point of 
this process is from the thinned rate Cp process combined with the rate 1 
process is {Cp + 1)/{C + 1). Therefore, the number of points generated in the 
rate C -F 1 process has a geometric distribution with mean {C + 1)/ {Cp + 1). 
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Each of these points has a C'/(C' + 1) chance of coming from the rate C 
process initially, and so requires a coin flip. Therefore, combining these effects 
gives an expected number of coin flips of [C/{C + 1)][(C + l)/{Cp + 1)] = 
C/{r + l). □ 


Now suppose that there is a known M < 1 /2 such th at r < M. Let 
BF(C) denote the Bernoulli Factory from Hubert ( to appe^ that flips a Cp 
coin using on average 9.5C(1 — M)~^ flips of the original coin. Consider the 
following algorithm. 


Small_r_lD_Bernoulli_Factory Input: C, M 

1) /3^1/{1-2M) 

2) Draw Y Logistic_Bernoulli_Factory(/5C') 

3) Draw H ^ Bern (1//9) 

4) liY = 0, then X ^ 0 

5) Elseif Y = 1 and B = 1, then X ■(— 1 

6) Else X ^ BF{I3C/{I3 - 1)) 

Lemma 3. Algorithm Small_r_lD_Bernoulli_Factory produces a Bernoulli 
distributed output with mean Cp < M < 1/2, and requires at most (on 
average) 

c I 

(1-2M)(1 + C'p) ^ ^ 

coin flips to do so. 

Note that for small p and M, this running time is to hrst order just C. 

Proof. First show correctness. Let Ai be the event that Y = 1 and H = 1 in 
the algorithm (in which case line 5 sets X to be 1), and A 2 be the event that 
Y = 1, H = 0, and a call to BF(C/5/(/5 — 1)) returns a 1 (in which case line 7 
sets X to be 1). These are disjoint events, and the output of the algorithm 
is X = l(Hi) + 1 (^ 2 ). Therefore, 


19C 


1 


1-2M + Cp 


P(X = 1) =P(Hi)+P(H 2 ). 


The value of Y is the call to Logistic_Bernoulli_Factory(/3C'), and 
so P(X = 1) = flCp/il + flCp). For the Bernoulli B, P(i? = 1) = 1/fl. 
Therefore P(H) = P(y = 1)P(H = 1) = Cp/{1 + flCp). 















The output of BF(C/9/(/9 — 1)) is 1 with probability equal to /3Cp, so 


P(^) = 


l + ^Cp 


(1 - 1 // 9 ) 


PCp 

P-1 


Therefore 


P(X = 1) = Cp- 


+ Cp- 


= Cp 


pcp 


PCp 
1 + pCp 


= Cp 


1 + pCp ^l + pcp 

as desired. 

Now for the running time. Lemma [2] gives a running time of pC/{l + Cp) 
for line 1. Line 5 is executed with probability {P — l)Cp/{l + PCp). By the 
way 0 was chosen, CpP/{P — 1) < 1/2. Therefore, Theorem 1.1 of Hnberl 
fjto appear! ) gives that the call to BF(C'/5/(/5 —1)) requires at most 19[C/3/(/3 — 
1)] flips. Therefore, the total number of flips is on average at most 


C 


{1-2M){1 + Cp) 


+ Cp' 


19C 


1-2M + Cp 


□ 


This is close to Theorem [H but the constant of 19 in the second term is 
larger. To improve this algorithm, and eliminate the need for the call to the 
old BF algorithm, it is necessary to consider what happens for larger r. 


3 Large r algorithm 


In this section, the algorithm of the previous section is improved to allow for 
all r G [0,1 — e], wh ere e is arbitrarily close to 0. Along the way, the older 
9.5C'e“^ algorithm of Hubert ( to appear ) is improved to a 7.5Ce~^ algorithm. 

The hrst step is to build a random coin flip whose mean is slightly larger 
than r. If this coin is tails, return tails for r. If the coin returns heads, heads 
will be returned with probability close to 1. Otherwise, a new coin will need 
to be flipped. 


3.1 A coin flip with mean slightly larger than r 

Consider an asymmetric random walk on the integers O = {0,1,...,m}, 
where given the current state X^, the next state is either max{ 0 ,X 4 — 1}, or 
minjXt + l,m}. The transition probabilities are 

P(Xi+i = min{i + l,m}\Xt = i) = Pr, P(Xi+i = max{z - 1,0 }jX* = i) = 


9 

















where Pr + Qr = 1- This is also called the Gambler’s Ruin walk. 

The following facts about this well known process will be helpful. 


Fact 7. Suppose Pr 7 ^ qr and T = inf{f : Xt G Then 


P(Xr = m) = 
E[T] = 

Fact 8. Suppose Xq = m, 
E[T] < m/{qr -Pr)- 


1 - {qr/Pr)^° 

1 - {qrlPvY 

Xq m 

qr Pr qr Pr 

Pr < qr, and T 


■ P(X'r = m). 
= inf{f : Xt 


( 2 ) 

(3) 

0}. Then 


Consider the following Bernoulli factory that begins a Gambler’s Ruin 
walk starting at state 1 , and returns heads if the state reaches 0 before it 
reaches m. 


A Input: m, C 

1) s t— 1 

2) While s G {1, 2, ..., m — 1} 

3) B -k— Logistic_Bernoulli_Factory(G) 

4) s^s-2B + l 

5) Return l(s = 0) 

Lemma 4. The output of k is a Bernoulli with mean r(l — — r™). 

Proof. Logistic_Bernoulli_Factory outputs a Bernoulli that has mean 
r/(l + r). Hence Pr = 1/(1 + r), qr = r/(l + r), and qr/Pr = r. From Fact [71 
in line 5 that makes P(/ = 0) = 1 — (1 — r)/(l — r™) = (r — r™)/(I — r™) = 
r(l-r™-^)/(l-r”^). □ 

Lemma 5. The expected number of coin flips used by A is at most C{m — 1). 

Proof. As in the last proof pr = I/{)- + r) and qr = r/(l + r). So 

r 1 1 — r 

1 + r 1 + r 1 + r’ 

and (qr/Pr) = r. Using (|3]), 


E[T] 


1 + r 
1 — r 


m- 


1 — r 

_ ^.m 



(1 + r) 


m 


1 

1 — r 
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Let /(r) = m /(1 — r™) — 1/(1 — r). Then it holds that /(r) < m — 1 for 
all r G (0,1) and m > 1. Note that for all r G (0,1): 


(1 - r™) 

/(r) < m — 1 ^ m{l — r) -< (m — 1)(1 — r^) 


^ m — mr — 1 — r — • • • — r 


m —1 


< m — 1 — r'^{m — 1) 


^ {m — l)r™' <r + r"^ + -- - + r'" ^ + mr 


<J=>m-l<r^ - 

Since r G (0,1), r*"™' > 1 for all i G {1,..., m — 1}, so the right hand side is 
strictly greater than the left hand side. Note /(O) = m —1, so for r G [0,1 —e], 
the function is at most m — 1 . 

Each call to line 3 requires on average (7/(1 + r) time by Lemma [2J 
Together, the overall number of steps (on average) is at most 


(1 + r)(m 


1 ) 


C 


1 + r 


C{m — 1). 


□ 


3.2 After the flip 


Here is how A can be useful. Using A, it is possible to generate a Bernoulli 
random variable that is 1 with probability 


P/3 


1 - 

' 1 - (/3r)’" 


for any constant /9 > 1. By choosing f5 large enough, pp > r. Note that 
P 13//3 <r. So r G [p/ 3 // 3 ,P/ 3 ]. 

So the algorithm works as follows. First flip a p/ 3 -coin. If it is tails, then 
return tails for the r-coin as well. If it is heads, then flip a (l//9)-coin. If 
that is heads as well, return heads for the r-coin. Otherwise, flip a p'-coin, 
and return the same value for the r-coin. 

For this algorithm to work, p' must satisfy: 


r = ^ + p'pf^{l - ( 1 //?)). 

Solving for p' gives 


! 

p = 


1 


(^r)™-i 

1 -I- (/5r)^ -1--h (/3r)"^“2 


1 


so the next step of the algorithm is figuring out how to generate a p'-coin. 
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3.3 Generating a p'-coin 


Fortunately, we do not have to actually generate a p'-coin for all possible 
values of (3, as we are allowed to choose the value of [3 to use, as long as 
> r for our choice of (3. Let 


so {(3 — 1) ^ = m — 1. Then 

/ _ jm - l){f3r)^-^ 

^ l + (/3r) + --- + (/?r)—2- 

Note that p' < 1 which gives that pp>r for this choice of (3. 

The algorithm for generating a p'-coin for p' as in (|1]) will be called B here, 
and is shown graphically in Figure [T] Notice that if the first flip is heads and 
the second flip is tails, then our problem has changed to the same problem, 
but with m reduced to m — 1. 



Figure 1: A graphical illustration of Algorithm B. 

To utilize this procedure, it is necessary to be able to generate a {(3r)'^~‘^/ (1+ 
■ ■ ■ + (/ 9 r )™'“2 coin. Fortunately, this can be accomplished fairly qnickly using 
the Gambler’s ruin chain from earlier. 
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High_Power_LogiStic_BF 

Input: Output: X ~ Bern((/9r)™/(l +-h (/Sr)™)) 

1) s ^ 1 

2) While s e {1,..., m} 

3) Draw B Logistic_Bernoulli_Factory(/9C') 

4) s^s + 2B-l 

5) X l(s = m + 1) 


Lemma 6. The output o/High_Power_Logistic_BF has distribution 6ern((/9r)™/(l+ 
■ • ■ + (/Sr)™)). The expected number of coin flips used is at most (3C/{I — (3r). 

Proof. This is a Gambler’s ruin where p = flr/{l + flr) and g = 1/(1 + flr), 
so q/p = l/(/9r). Hence from Fact [71 


(s = m + 1 ) = 


l — {l/{flr)fl (/3r)™-(l —/3r) 


{flr) 


1 - (l/(/5r))"^+i l-(/3r)"^+i IH-h (/5r)' 


Also from Fact [71 note q — p = (1 + flr)/{l — flr), so if T is the number 
of times line 3 is called, 


E[T] = 


1 + flr 
1 — /Sr 


1 — (m + 1) 


(/Sr)”^(l — flr) 
1 - (/5r)™+i 


< 


l + flr 
1 — flr' 


Each call to Logistic_Bernoulli_Factory takes time flC/{l + flr), so 
the overall number of coin flips (on average) is at most flC/{l — flr). □ 


In pseudocode, algorithm B looks like this. 


B Input: e,m,/3,C Output: X 

1) X ^ 0.5 

2) While X ^ {0,1} 

3) Draw Bi ^ Linear_Bernoulli_Factory(l — (1 — e)fl,fl ■ C) 

4) If Hi = 0 then X ^ 0 

5) Else 

6 ) B 2 High_Power_Logistic_BF(m — 2, fl, C) 

7) If H 2 = 1 then X ^ 1 

8 ) Else m •(— m — 1 
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The most important thing to note here is that like many perfect simn- 
lation algorithms, this method employs recnrsion. We do not yet have an 
algorithm for completing line 3! However, this algorithm B can be nsed as 
a snbrontine to create snch an algorithm, and then this snbrontine will call 
the hnished algorithm. 

Lemma 7. The output ofB has distribution 

Bern{{m - l)(/3r)™-V(l + ■ • • + (/5r)™-")). 


Proof. The proof is by indnction. When m = 2, if i?i = 1 then B 2 ~ Bern(l), 
so X = 1 with probability (dr as desired. 

Now snppose that the resnlt holds for m, consider m + 1. Then 


P(X = 1)= (3r 
(m 


{(dr\ 


im— 2 


_H-h (/3r)™-2 


IH-h i/dr) 


m-2 ’ 


1 + ■ • • + {(dr)^-^ 
iH-h (/5r)™-2 


(m — 2)(/3r)™' ^ 
iH-h (/5r)”"-3 


completing the indnction. 


□ 


3.4 The new linear Bernoulli Factory 

With these preliminaries in place, the overall algorithm is as follows. 


Linear_Bernoulli_Factory Input: e,C Output: B 

1) 

m [4.56 ^] + 1, 1 + l/(m — 1) 

2) 

Bi A(m, (d ■ C) 

3) 

If Hi = 1 

4) 

Draw B 2 Bern(l//3) 

5) 

If H 2 = 1 then H 1 

6) 

Else 

7) 

Draw B ^ B(m, (d, C) 

8) 

Else H 0 


Now consider the expected nnmber of coin flips nsed by the algorithm. 
As will become clear in the proofs of Lemmas [H] and (TD] below, making m = 
0(e“^) is the correct choice. That leaves the choice of constant np to us, and 
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the constant of 4.5 from Line 1 was chosen to make the overall running time 
as small as possible. 

This algorithm calls A and B. Line 2 of B needs to draw a Bern (/dr) random 
variable. The best way to draw these random variables is to recursively call 
Linear_Bernoulli_Factory. In order to ensure that this back in forth calling 
eventually comes to a halt with probability 1, it is easiest to bound the total 
expected number of calls to Linear_Bernoulli_Factory. 

Lemma 8. The expected number of calls to Linear_Bernoulli_Factory is 
at most l.f. 

Proof. Let mi be the value of m in the first call to Linear_Bernoulli_Factory, 
and /di = 1 + l/(mi — 1). From this hrst call there is a chance of calling B, 
which in turn calls Linear_Bernoulli_Factory with m 2 and /d 2 - Each of 
those second generation calls might call a third generation, and so on. To 
bound the expected number of calls to Linear_Bernoulli_Factory sum over 
all possible calls of the probability that that call is executed. Let iVj denote 
the number of ith generation calls. 

The expected number of calls in the hrst generation is 1. Consider a call 
in the second generation. In order for that call to be made, there must have 
been a call to B from the hrst generation, and all prior second generation calls 
from line 3 of B must have had Bi = 0. The number of times the while loop 
in B is executed is stochastically dominated by a geometric random variable 
with mean 1/(1 — /Sir). Since 

1-/Sir > l-(l + l/[4.5e-il)(l-e) (5) 

= (7/9)e+(2/9)e^ (6) 

the number of calls made is bounded (in expectation) by (9/7)e“^. 

But before B is even called, hrst it must have held that Bi = 1 and B 2 = 0 
in lines 2 and 4 of the hrst generation call to Linear_Bernoulli_Factory. 

The probability that a call to B is made is at most 

(1 - l/ft)/3r(l - (Ar)"-‘)/(l - (I3,rn < 

m — 1 

So the expected number of calls to Linear_Bernoulli_Factory in the 
second generation is bounded by 

E[iV 2 |iVi] < iVi[l/(mi - l)](9/7)e-' < (2/7)iVi. 
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This step forms the basis of an induction that gives ]E[iVj] < (2/7)®iVi = 
(2/7)h Therefore < 1/(1 - 2/7) = 1.4. □ 

Lemma 9. The output B o/Linear_Bernoulli_Factory has B ~ Bern{r). 

Proof. Line 2 of B requires a draw Bi ^ Bern(/9r). Suppose that for the 
hrst L times this line is called, the Linear_Bernoulli_Factory is called to 
generate this random variable. Then, from the L + 1st time onwards, an 
oracle generates the random variable. 

We show by strong induction that Linear_Bernoulli_Factory generates 
from Bern(r) for any hnite M. The base case when M = 0 operates as follows. 
Lemma [7] immediately gives in this case that a call to B returns a random 
variable with distribution Bern((m—l)(/5r)”^“^/(1 + - • •+(/9r)™'“^)). LemmalU 
gives that Bi from line 2 has distribution Bern((/9r)™'/(l + ■■■ + 

Putting this together gives 


P(5 = 1) 


= {/3r) 
= Wr) 
= r 


m—1 


1 - 

1 - (/dr)"* 

1 - (/dr)”*-^ 
1 - (/dr)"* 
(1 - (/dr)”*-^) 
1 - (/dr)"* 


1 / 1\ (m —l)(/dr)"* ^ 

1 + ■ ■ ■ + (/3r)”*-2 
1 , 1 (/dr)"*-^ 

^ /d IH-h (/dr)”*-2_ 

1 + ■ ■ ■ + (/dr)"*-^ 

IH-h 


= r. 


This is the rare induction proof where the base case is just as hard as 
the induction step. Suppose it holds for L, and consider what happens for 
call limit L + 1. Then the hrst call to Linear_Bernoulli_Factory might call 
B, which might call Linear_Bernoulli_Factory. But the hrst such call has 
used up one call, so only has L +1 — 1 calls remaining, so by strong induction 
each returns the correct distribution. Hence Lemma [7] holds, and the hrst 
call returns the correct distribution by the same argument as the base case. 

Let N be the random number of calls to Linear_Bernoulli_Factory 
needed by the algorithm. Then let B be the output when N is unbounded, 
and Bm be the output when a limit on calls equal to L is in place. Then 

P(H = 1) = P(5 = l,N <L)+¥{B = l,N > L) 

= P(Hi = l,N <M) + P{B = l,N > L) 

= ¥{Bl = 1) - P{Bl = 1, at > L) + P(H = 1, iV > L). 
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Both P(-Bl, N > L) and P(-B = 1, > L) are bounded above by P(A^ > L). 

Since by the last lemma E[iV] < 1.4, limi^ooP(iV > L) = 0. The only way 
this can hold for all L is if F{B = 1) = P(i?L = 1) for all L, so B has the 
correct distribution. □ 

Lemma 10. Linear_Bernoulli_Factory uses on average at most 7.67Ce~^ 
coin flips to generate B ~ Bern{r). 

Proof. From the proof of Lemma [HI the expected number of calls to the fth 
generation of Linear_Bernoulli_Factory is bounded above by (2/7)*. 

From ([6]), at each successive generation of calls, e is being multiplied by a 
factor of at least 7/9. So an ith generation call to Linear_Bernoulli_Factory 
has an m value of at most |'4.5(9/7)*e“^] + 1, where e was the input for the 
0th generation. 

Coin flips occur during the call to A, and Lemma 0] bounds the ex¬ 
pected number of coin flips by C[4.5(9/7)*e“^]. So the expected total flips 
coming from the ith generation of Linear_Bernoulli_Factory is at most 
C'[4.5(18/49)*e-i + (2/7)*]. 

Now look at the flips coming from an ith generation call to B. This genera¬ 
tion is only called from an ith generation call to Linear_Bernoulli_Factory, 
of which there the expected number is at most (2/7)*. The call to B occurs 
with probability at most [fl — l)r, so at most r/m. The while loop inside 
B is run (on average) at most e~^ times, each of which could make a call 
to High_Power_Logistic_BF. By Lemma [6] this requires at most flCe~^ coin 
flips. So the total number of coin flips from an ith generation call to B is at 
most 


(2/7)*[r/m][/3C'e“^(9/7)*+^]e"^ < (9/7)4.5“^(18/49)*Ce"^ 

Summing over these flips and the ones from Linear_Bernoulli_Factory 
gives a total sum of 

(469/62)C'e"^ < 7.57C'e"^ 

coin flips on average. □ 

3.5 Small r 

Now that a recursive algorithm for large r has been built, a recursive analogue 
for Small_r_lD_Bernoulli_Factory works as follows. 
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Small_r_Bernoulli_Factory Input: C, M 

1) /3^1/(1-2M) 

2) Draw Y Logistic_Bernoulli_Factory(/9C') 

3) Draw-B ^ Bern (1//9) 

4) If y = 0, then X ^ 0 

5) Elseif y = 1 and i? = 1, then X 1 

6) Else X ■(— Linear_Bernoulli_Factory(C'/5/(/3 — 1)) 

Lemma 11. Algorithm Small_r_Bernoulli_Factory produces a Bernoulli 
distributed output with mean Cp < M < 1/2, and requires at most (on 
average) 

c I 

(1-2M)(1 + C'p) ^ ^ 

coin flips to do so. 

Proof. The proof is essentially the same as that of Lemma |3l □ 


C 


15.2 


1-2M + Cp 


4 Lower bound 


To see why it is unlikely that a method that uses fewer than Ce~^ coin flips 
can be constructed, consider building an unbiased estimate of p. 

The standard estimate is to generate Xi,... ,X„ iid Bern(p), and then 
use the sample average Pn = (Xi + • • • + Xflj/n as an unbiased estimate of 
p. This estimate is unbiased, and has variance p(l —p)/n. 

Now consider the estimate Y/C, where Y ~ Bern(C'p). Then E[y/C'] = 
Cp/C = p so this estimate is also unbiased, and the variance is Cp{l — 
Cp)/C‘^ = p{l — Cp)/C. Therefore, this estimate that used one draw from 
Bern(C'p) has the variance of the estimate that used n = C{1 — p)/{l — Cp) 
draws from the p-coin. 

The Cramer-Rao lower bound (see, for instance Bickel and Doksum (1973)) 
on the variance of an unbiased estimate of p is 

P (1 -P) 
n 

That is, any unbiased estimate that uses up to n flips of the p-coin must have 
variance at least p(l — p)/n. That immediately gives that any algorithm for 
generating a Cp-coin that uses a deterministic number n of coin flips must 
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have n > C{1 — p)/{I — Cp). Of course, this does not quite apply to a 
Bernoulli Factory, because here a random number of coin flips is used. 

However, it is strong evidence that 0(1 — p)/(l — Cp) is a lower bound 
on the expected number of coin flips needed by an algorithm. 


5 Multivariate Bernoulli Factory 

This new algorithm was designed for the single coin problem, in this section 
consider generating a coin flip whose probability of heads is the sum of the 
probability of heads on two different coins each of which has an unknown 
probability of heads. Unlike the single coin flip, there is no immediate ap¬ 
plication, however, it is useful to know that the single coin algorithm can be 
easily generalized to solve this problem should the need arise. 

More generally, the goal is now to generate a coin flip with probability 


r = Cipi H-h CkPk 


of heads, where r is bounded away from 1 , using as few flips of the coins 
as possible. When k = 1, this is the linear Bernoulli factory studied in the 
previous sections. Formally, a multivariate Bernoulli factory is dehned as 
follows. 

Definition 3. Given a computable function / : [0,1]^ —)■ [0,1], a multivariate 
Bernoulli faetory is a computable function 

.4:|0,l]x ({ 0 , 1 }W I)‘^{ 0 , 1 }, 

such that if U ~ Unif([0,1]) and the Xij are independent random variables 
with (Vi G {!,..., k}){'ij G {1, 2 ,.. .})(Xjj ~ Xi), then the following prop¬ 
erties hold. 

1. There exist random variables (Ti,..., T^) G {1, 2,.. .}^ such that the 

value of Vl(U, ..., only depends on the values of 

• • •) all (ti,..., 4 ), the event (Ti,..., T^) = 

(ti,..., tk) is measurable with respect to ..., {Xk^ijl'L^. 

2. A{U, {Xi^ijZi, • • •, ~ Bem(/(pi,p2, • • ■,Pk))- 

Call Ti + ■ ■ ■ + Tk the running time of the algorithm. 
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The key to the single coin algorithm was generating a random variable 
that was exponential with rate parameter r. In the single coin case, this was 
done by generating a Poisson process of rate C, then thinning. 

For the multivariate coin case, consider generating k independent Poisson 
point processes Pi,..., P^, where P* has rate Ci. Then thin each process P* 
with coin i to obtain a Poisson point process of rate CiPi. The union of these is 
a new Poisson point process of rate r, and the rest of the algorithm operates 
as before. The proofs of all the lemmas and Theorem |2] then proceeds in 
exactly the same way as given earlier. 
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